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This paper presents an on-the-fly uniformization technique for the analysis of time-inhomogeneous 
Markov population models. This technique is applicable to models with infinite state spaces and 
unbounded rates, which are, for instance, encountered in the realm of biochemical reaction networks. 
To deal with the infinite state space, we dynamically maintain a finite subset of the states where most 
of the probability mass is located. This approach yields an under-approximation of the original, 
infinite system. We present experimental results to show the applicability of our technique. 

1 Introduction 

Markov population models (MPMs) are continuous-time Markov processes, where the state of the sys- 
tem is a vector of natural numbers (i.e., the populations). Such models are used in various application 
domains: biology, where the state variables describe the population sizes of different organisms, queue- 
ing theory, where we model a state as a vector of queue occupancies, chemistry, where the state variables 
represent the amount of molecules of different chemical species, etc ifTTIl . 

Besides the expectations and variances of the different populations, the probabilities of certain events 
occurring can be of interest when studying MPMs. It may be necessary to know the probability of the 
extinction of a species, the probability that a population reaches a certain threshold, or even the full 
distribution of the MPM at a certain time-point, for instance to calibrate model parameters. 

Many Markov population models have infinitely many states. In the case of biological or chemical 
applications, we normally cannot provide hard upper bounds for population numbers and in the field of 
queueing theory it may be interesting to consider unbounded queues. The evaluation of infinite MPMs 
through numerical H or statistical Q analysis has been well-studied for time-homogeneous models 
where the dynamics of the system are independent of time. In [4J the state space of the model is generated 
and truncated on-the-fly during the transient solution, that is, during a certain time interval only states 
that are relevant at that time are considered. Thus, states are added at a certain step and dropped at a 
later time when they become irrelevant. A similar technique is proposed in |[3l for the solution of time- 
homogeneous discrete-time Markov chains. Note that this is different from on-the-fly techniques for the 
computation of steady-state probabilities where the relevent part of the state space is generated but states 
are never dropped as time progresses llT4l . 

Many Markov models are time-inhomogeneous, that is, their dynamics change over time. For in- 
stance, when modeling an epidemic, we may have to take into account that infection rates vary season- 
ally. For traffic models, time-dependent arrival rates can be used to model the morning and evening rush 
hours. In cellular biology we see that reaction propensities depend on the cell volume, which waxes and 
wanes as the cell grows and divides. The class of finite time-inhomogeneous Markov models has also 
been studied in recent years ll2ll5l [T3i . 
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In this paper, we develop a numerical algorithm to approximate transient probability distributions 
(i.e., the probability to be in a certain state at a certain time) for infinite time-inhomogeneous MPMs. We 
consider MPMs with state-dependent rates and do not require the existence of an upper-bound for the 
transition rates in the MPM. 

Our algorithm is based on the uniformization technique, which is a well-known method to approxi- 
mate the transient probability distribution of finite time-homogeneous Mai^kov models |[T0l l9l. Recently, 
two adaptations of uniformization have been developed. These adaptations respectively approximate 
the transient probabilities for finite time-inhomogeneous [2] and infinite time-homogeneous HI Markov 
models. Our algorithm combines and refines these two techniques such that infinite time-inhomogeneous 
MPMs with unbounded rates can be tackled. We present two case studies to investigate the effectiveness 
of our approach. 

2 Markov Population Models 

Markov chains with large or even infinite state spaces are usually described by some high-level model- 
ing formalism that allows the generation of a (possibly infinite) set of states and transitions. Here, we 
use transition classes to specify a Markov population model, that is, a continuous-time Mai^kov chain 
(CTMC) {X{t),t > 0} with state space 5 = Z'| = {0, 1, . . .}", where the i-th state variable represents the 
number of instances of the i-th species. Depending on the application area, "species" stands for types of 
system components, molecules, customers, etc. The application areas that we have in mind are chemical 
reaction networks, performance evaluation of computer systems, logistics, epidemics, etc ifTTI . 

Definition 1 (Transition Class) A transition class x is a triple {G,w,a) where G C Z'^ is the guard, 
w € Z" is the change vector, and a . G x M>o — ?• M>o is the time-dependent rate function. Moreover, for 
any x € we have that x ^G implies x-\-w ^ Z']_. 

The guard is the set of states where an instance of T is possible, and if the cuiTcnt state is x E G then 
;c + w € is the state after an instance of T has occurred. The rate a{x,t) determines the time-dependent 
transition probabilities for an infinitesimal time-step dt 

Pr{X{t + dt) =x + w \X{t) =x) = a{x,t)-dt + o{dt), 

where o is a function such that o(0) = and limi,^oo{h)/h = 0. 

A CTMC X can be specified by a set of m transition classes Ti , . . . , as follows. For j € { 1 , . . . , m}, 
let Tj = {Gj,Wj, aj). For each t € M>o we define the generator matrix Q{t) of X such that the row that 
describes the transitions of a state x has entry aj{x,t) at position Q{t)x.x+w whenever x G Gj and zero 
otherwise. Moreover, the diagonal entries of 2(0 are the negative sums of the off-diagonal row entries 
because the row sums of a generator matrix aie zero. We assume that each change vector wj has at 
least one non-zero entry. To simplify the presentation we assume that all change vectors are distinct. 
We remark that X is called time-homogeneous when Q{t) is equal for all t. Otherwise, X is called time- 
inhomogeneous. 

Example 1 We consider a simple gene expression model for E. coli cells l[I6\l . It consists of the tran- 
scription of a gene into messenger RNA (mRNA) and subsequent translation of the latter into proteins. A 
state of the system is uniquely determined by the number of mRNA and protein molecules, that is, a state 
is a pair {xm,xp) G Z^. We assume that initially there are no mRNA molecules and no proteins in the 
system, i.e., Pr{X{0) = (0,0)) = 1. Four types of reactions occur in the system. Let j G {1, ... ,4} and 
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Tj = {Gj,Wj,(Xj) be the transition class that describes the j-th reaction type. We first define the guard 
sets Gi , . . . , G4 and the change vectors wi , . . . , W4. 

• Transition class Zi models gene transcription. The corresponding stoichiometric equation is — ?> 
mRNA. If a Zi-transition occurs, the number of mRNA molecules increases by one. Thus, w\ = 
(1,0). This transition class is possible in all states, i.e., G\ = Z^. 

• We represent the translation of mRNA into protein by T2 (mRNA — )■ niRNA+P). A Z2-transition is 
only possible if there is at least one mRNA molecule in the system. We set G2 = {{xm,xp) € | 

> 0} and wj = (0, 1). Note that in this case mRNA is a reactant that is not consumed. 

• Both mRNA and protein molecules can degrade, which is modelled by T3 and T4 ( mRNA — )> and 
P — )■ 0j. Hence, G3 = G2, G4 = {{xm,xp) e | xp > 0}, W3 = (—1,0), andw4 = (0,-1). 

Let k\,k2,k-i,k4 be real-valued positive constants. We assume that transcription happens at rate 
ai{xM,xp,t) = k\ ■V{t), that is, the rate is proportional to the cell volume V{t) M 71/ . The (time- 
independent) translation rate depends linearly on the number of mRNA molecules. Therefore, Ot2{xM,xp,t) 
= k2 -xm. Finally, for degradation, we set a3{xM,xp,t) = k^-XM and a4{xM,xp,t) = k^ -xp. 

We now discuss the transient probability distribution of a MPM. Let S be the state space of X and let 
the transition function P(f,f + A) be such that the entry for the pair {x,y) of states equals 

P{t,t + ^)^^=Pr{X{t + ^)=y\X{t)=x), ?,A>0. 

If the initial probabilities Pr(X(0) =x) are specified for each x G 5, the transient state probabilities 

p^'\x) := Pr(X{t) = x), are given by 

We assume that a transition class description uniquely specifies a CTMC and rule out "pathological 
cases" by assuming that the sample paths X{t) are right-continuous step functions. In this case the 
transition functions are the unique solution of the Kolmogorov backward and forward equations 

jPM = Q{t)-P{to,t) (1) 
jPM=PM-Q{t), (2) 

where <to<t. Multiplication of Eq. with the row vector p^'") ^^jth entries p^'^^x) gives 

^.pi^)=pi^).Q(t). (3) 

If S is finite, algorithms for the computation of p^'^ are usually based on the numerical integration of the 
linear system of differential equations in Eq. ^ with initial condition p^'^\ Here, we focus on another ap- 
proach called uniformization that is widely used for time-homogeneous Markov chains ifTQl. It has been 
adapted for time-inhomogeneous Markov chains by Van Dijk Q and subsequently improved |[T3l l2l. 
The main advantage of solution techniques based on uniformization is that they provide an underapprox- 
imation of the vector p'^'^ and, thus, provide tight error bounds. Moreover, they are numerically stable 
and often superior to numerical integration methods in terms of running times lITSl . 
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3 Uniformization 

Uniformization is based on the idea to construct, for a CTMC X, a Poisson process N{t),t >0 and a 
subordinated discrete-time Markov chain (DTMC) Y (/) , / G N such that for all x and for all t 

Pr{X{t)=x)=Pr{Y{N{t))=x). (4) 

Since Poisson process and DTMC Y are independent, the equation above can be written as 

oo 

Pr{Y{N{t))=x) = £Pr(F(/) =x)Pr{N{t) = i) . (5) 

i=0 

For a finite time-homogeneous MPM with state space S the rate A of the Poisson process (also called 
the uniformization rate) is chosen to be greater than or equal to the maximal exit-rate appearing in X 



A > max Y\ aj (x) . 
For the DTMC Y we find transition probabilities 

OC " (x) 

Pr{Y{i+\) =x+Wj I Y{i) =x) = 

When X is time-inhomogeneous, Arns et al. Ill suggest to define the time-dependent uniformization rate 
A{t) of the inltomogeneous Poisson process (IPP) A'^ as 

m 

A{t) > max ^ aj{x,t). (6) 

■'■^^ y=i 

(X '(x 

For the (time-dependent) transition probabilities of the DTMC Y we then have that is the proba- 

bility to enter state x + wj from state ;c if a state-change occurs at time t. Arns et al. prove that Eq. (|4]) 
holds if the aj are (right or left) continuous functions in t and if S is finite (see Theorem 7 in lH). Here, 
we relax the latter condition and allow S to be infinite. If sup^^^^y aj{x,t) < o° during the time interval 
of interest, the proof of Eq. ^ may be expected to proceed along similar lines. In our case, however, 
sup^g^^y o:y(A:,f) = oo and then the Poisson process A'^ is not well-defined as its rate must be infinite 
according to Eq. Therefore, the infinite state space has to be truncated in an appropriate way. 



3.1 State Space Truncation 

We consider a time interval [t,t + A) of length A, where the transient distribution at time t, p^'\ of 
the infinite time-inhomogeneous MPM X is known. We now wish to approximate the transient distri- 
bution at time t + A, p^'^^\ We assume that p^'^ has finite support St,o. Define Pr{N{t,t + A) = i) = 
Pr{N{t + A) — N{t) = i) as the probability that A^ performs / steps within [t,t + A). For a fixed positive 
e <S 1, let R and the rate function A be such that St,R is the set of states that are reachable, with probability 
greater than or equal to 1 — £, from the set St a in the time-interval [t,t + A) within at most R transitions, 
i.e. 

R 

Y,Pr{N{t,t + A) = i)>l-e. (7) 

(=0 
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Furthermore, we have that the rate of A'^ at time t' G [t,t + A) must satisfy 

m 

A(?') > max V ay(x,f')- (8) 

Note that A(?') is adaptive and depends on t', t. A, Stfi, and R as opposed to Arns et al. where A(f') 
depends only on f', t, and A because they consider finite state spaces. 

Finding appropriate values for A and R is non-trivial as A(?') determines the speed of the Poisson 
process A'^ and thereby influences the value of 7?. On the other hand, R determines the size of the set St^R 
and thus influences A(f'). We discuss how to find appropriate choices for A and R given the set St^o in 
Section O] 

Assume that we find A and R with the above mentioned properties and define A(f') as in Eq. ([8]l. 
Then, for all x G 5, we get an £-approximation 

R 

Pr{X{t+A) =x) > £ Pr{Y{i) =xAN{t,t+A) = i) , (9) 

(=0 

where Y has initial distribution p^'\ The probabilities Pr{Y{i) = x AN{t,t + A) = i) can now be approx- 
imated in the same way as for the finite case [21. 

From Eq. ^ we see that it is beneficial if R is small, since this means fewer probabilities have to be 
computed in the right-hand side of Eq. Note that the truncation-point R is small when the uniformiza- 
tion rates A(?') are small during [t,t + A) because if A'^ jumps at a slower rate then Pr{N{t,t + A) > /) 
becomes smaller. Thus, it is beneficial to choose A{t') as small as possible while still satisfying Eq. ([8]). 



3.2 Bounding approach 

Let p^'^^\x) denote the right hand side of Eq. Q, i.e., the approximation of the transient probability of 
state X at time t + A. We compute this approximation with the uniformization method as follows. The 
processes Y and are independent which implies that 

Pr{Y{i) =xAN{t,t+A) = i) = Pr{Y{i) =x) ■ Pr{N{t,t+A) = i) . 

The probabilities Pr{N{t,t + A) = /) follow a Poisson distribution with parameter A{t,t + A) • A, where 

A(f,f + A) = l/;+^A(0^?'. 

For the distribution Pr{Y{i) =x), Arns et al. suggest an underapproximation that relies on the fact that 
for any time-point G [f,? + A) we have: 

Thus, for / G {1,2, . . . ,/?}, we iteratively approximate Pr{Y{i) =y) as 

Pr{Y{i)=y)> £ Pr{Y{i-\)=x)-Uj{x,t,t+A)+Pr{Y{i-\)=y)-uo{y,t,t+A). (10) 

x,j:y=x+Wj 

Here, x ranges over all direct predecessors of y and the self-loop probability uo{y,t,t + A) of 3^ is given 
by 

uo{y,t,t + A)= mm [l- L ^ ) . 

t'elt,t+A] \ j=i I I 
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Note that often we can split aj{x,t') into two factors Ay(f') and rj{x) such that aj{x,t') = Xj{t') ■ rj{x) 
for all t'J, Thus, the functions Xj : M>o — >• IR>o contain the time-dependent part (but are state- 
independent) and the functions rj : S ^ M>o contain the state-dependent part (but are time-independent). 
Then each minimum defined above can be computed for all states by considering 



mm 



In particular, if Ay and A are monotone, the above minimum is easily found analytically. 

The approximation in Eq. (ITOl ) implies that for the time interval [t,t +A), we compute a sequence 
of substochastic vectors v^^\v^'^\ . . . ,v^^'> to approximate the probabilities Pr{Y{i) =y). Initially we 
start the DTMC Y with the approximation p^'^ =: v^''' of the previous step. Then we compute v^'^^^ 
from v^'' based on the transition probabilities Uj{x,t,t + A) for / e {0, 1, . . . ,R}. Since these transition 
probabilities may sum up to less than one, the resulting vector v^'+^' may also sum up to less than one. 
Since, for the computation of p'^^, we weight these vectors with the Poisson probabilities and add them 
up the underapproximation p'^^ contains an additional approximation error. In general, the larger the 
time-period A, the worse the underapproximations Uj{x,t,t + A) are and thus the underapproximation 
p'~^^ becomes worse as well. We illustrate this effect by applying the bounding approach to our running 
example. 

Example 2 In the gene expression of Example |7] the time-dependence is due to the volume and only 
affects the rate function ai of the first transition class. The time until an E. coli cell divides varies 
widely from about 20 minutes to many hours and depends on growth conditions. Here, we assume a cell 
cycle time of one hour and a linear growth Thus, if at time t = we consider a cell immediately 
after division then the cell volume doubles after 3600 sec. Assume that A < 3600. Then, ai{x,t') = 
k\- (1 + jgQQ ) for all X € 5. Assume we have a right truncation point R such that 

t' 

A{t') = maxk\ • (1 + — — ) + (^2 + ^3) • -^r + ^4 • xp 
XR,xp 3600 

where xr and xp range over all states {xp,xp) € Sq^r and Eq. (|7]l holds. Then we find, for each time-point 
t' G [0,A), the same state for which the exit-rate (Xo{x,t') := Y!J=i OCjix,t') is maximal, since the only time- 
dependent propensity is independent of the state-variables. Let (x^^^jX™^") denote this state. In general 
this is not the case, for instance in the realm of chemical reaction systems we have that the propensities 
of bimolecular reactions (reactions of the from A+S — > ...) are dependent both on cell-volume and 
the population numbers. For such a system we may find that different states have the maximal exit-rate 
within the time-frame [0,A). We discuss how to overcome this difficulty in Subsection \4.2\ The transition 
probabilities of the DTMC Y are now defined as 

. ai(xR,xp,?') ai(x,0) k[ 
Mi(xR,xp,0, A) = mm — — 



t'em A{t') A(0) k[ + {k2 + k^)-xf^^ + k^-xf' 

and, for 7 G {2,3}, 

. aj{xR,xp,t') . kj-XR kj-XR 

Uj{xR,xp,0,A) = mm — — — = mm — 



,'em A{t') r'e[o,A]A(A) k[ ■ (l + ^) + {k2 + k^) ■ x'^'^ + k^ ■ x'^'^ ' 



'Note that this decomposition is always possible for chemical reaction networks where the time-dependence stems from 
fluctuations in reaction volume or temperature. 
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Figure 1: Illustration of the state space truncation approach for the two-dimensional case. Given the 
distribution p(') with support Stfl, a truncation point R and a time-step A, we compute in the first step the 
distribution p^'+^^ with support St.R = ^f+A.o- For the next step we consider the set St+A,R- 



k4 -Xp 



Ua(xr,Xp,0,A) = 7 — 

^'r (1 + 3^) + (^^2 + ^3) + ^4 •X^^'' 

For the self-loop probability we find: 

, n A\ -1^4^ aj{xR,xp,t')\ ( ^ aj{xR,xp,t') 
uo(xR,xp,0,A) = mm 1- > . , ,, = 1" max > , , ,, 

^ aj{xR,xp,A) ^ k\-{l + j^) + {k2 + k-i)-XR + k4-xp 



h A(A) " k\-i\ 



j=i A(A) k\-{\ + ^) + {k2 + k^)-xf-- + k4-xf--' 

We now calculate the fraction of probability lost during the computation of v^'^^^ from v('), i.e., 



1 ~ H My(xR,Xp,0,A) 



Po'' ' k[.{\ + ^) + {k2 + k^).^^^ + k,.xT^ k[ + ik2+k,).xf^^ + k,.xf^' 

{k2 +k3)- X^^" + k4 ■ X]?^" {k2 +k3)- X^''" + k4 ■ x"?^" 



max 



k[ + ik2+k3) -X-- + ^4 -^r" K-i'^ + mo) + (^2 + /^3) ■Xf'^ + k^-x] 

For A = we have a probability loss of and for A > we can see that the probability loss increases 
with increasing A. 



3.3 Time-stepping approach 

Given that a large time horizon may lead to decreased accuracy, Arns et al. lH suggest to partition the 
time period of interest [0,?niax) in steps of length A. In each step, an approximation of the transient 
distribution at the current time instant, p^'\ is computed and used as initial condition for the next step. 
The number of states that we consider, that is, |5, ^| grows in each step. The probabilities of all remaining 
states of S are approximated as zero. Thus, each step yields a vector ^('+^) with positive entries for all 
states X € St^R that approximate Pr{X{t + A) = x). The vector p^*+^^ with support St^p = St+Afi is then 
used as the initial distribution to approximate the vector pi'+^+^\ See Figure [T] for a sketch of the state 
truncation approach. Note that the chosen time-period A may vary for different steps of the approach. 
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It is easy to see that the total error is the sum of the errors in each step, where the error of a single step 
equals the amount of probability mass that "got lost" due to the underapproximation. More precisely, we 
have two sources of error, namely the error due to the truncation of the infinite sum in Eq. ^ and the 
error due to the bounding approach that relies on Eq. (flOl ). 

In lIU, Arns et al. give exact formulas for the first three terms of the sum in Eq. © (for / = 0, 1,2). 
Thus, if the approximation p^'^ of p^'^ is exact, then is an underapproximation due to the remaining 

terms in Eq. This implies that the smaller R becomes, the closer the eiTor will be to the eiTor 
bound £. On the other hand, a small truncation point means that only a small time step A is possible (see 
Eq. (|71)), which means that many steps are necessary until the final time instant ?max is reached. In order 
to explore the trade-off between running time and accuracy, we run experiments with different values 
for the predefined truncation point R that determines the step size A. We report on these experiments in 
Section |5] 

4 On-the-fly Algorithm 

As we can see in Figure [B the number of states that are considered to compute from p^'^ grows in 

each step, since all states within a radius of R transitions from a state in the previous set St a are added. 
This makes the approach infeasible for Markov models with a large or even infinite state space because 
the memory requirements are too large. Therefore, we suggest to use a similar- strategy as described in 
previous work i4j| to keep the memory requirements low and achieve faster running times. 

The underlying principle of this approach is to dynamically maintain a snapshot of the part of the 
state space where most of the transient probability distribution is located. We achieve this by adding 
and removing states in an on-the-fly fashion. The decision which states to add and which states to 
remove depends on a small probability threshold 5 > 0. After the computation of the vector v('+'' 
based on v^'-*, we set all entries in v^'+'^ to zero that have a probability less than 5. This significantly 
reduces the computational complexity since only parts of the transition probability matrix of Y have to 
be generated p] (for instance, we explore 360000 states at time instant t = 600 for the gene expression 
system of Example[T]if 5 = but only 5700 states are stored when 5 = 10^'^). Let 

and, for / G {1, . . . ,R} let S^'^ be the set of states that we consider to compute v('+^^ from v^'^. We remark 
that this also decreases the speed of the Poisson process A'^ since the sets 5f,o and St,R are smaller and thus 
the maximum in Eq. ^ is now taken over fewer states. We illustrate this effect in Figure [2] This effect 
is particularly important if during an interval in certain parts of the state space the dynamics of 

the system is fast while it is slow in other parts where the latter contain the main part of the probability 
mass. On the other hand, the threshold 5 introduces another approximation enw which may become 
large if the time horizon of interest is long. Moreover, if p is a bound for the error introduced by the 
above strategy of neglecting certain states, we can reserve a portion of the probability loss p • ^ for the 
interval [t,t + A) and repeat the computation with a smaller threshold 5 if more than the allowed portion 
of probability was neglected. 

The approximation that we suggest above is again an underapproximation and since the approxima- 
tions suggested in the previous sections ai^e also underapproximations, we ai"e still able to compute the 
total error of the approximation p^'^ of p^''^ as 

i-£^W(x). (11) 
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Support at time t 



X2 



Truncation for the first step 
I and approx. support of p''^^ 




t 

X2 



Truncation for the second step 




Figure 2: Illustration of the on-the-fly algorithm for the two-dimensional case. Given the distribution 
p^'^ with support 5, 0, a truncation point R and a time-step A, we compute in the first step the distribution 
p('+^) with approximate support St+A.o C 5, ^. For the next step we consider the set Si+a^r. 

Clearly, t' > t implies that the error at time t' is higher than the error at time t. For our experimental 
results in Section|5]we choose 5 = 10~^^ and report on the total eiTor of the approximation at time fmax- 



4.1 Determining the step-size 

Given an error bound e > 0, a time-point t, for which the support of p^'^ is Sf,o, and a time-point fmax for 
which we wish to approximate the transient probability distribution, we now discuss how to find a time- 
step A such that Eqs. ^ and dSJl hold. Recall that the probabihties Pr{N{t,t + A) = i) follow a Poisson 
distribution with parameter A{t,t + A) A, which we denote by (Ir a to emphasize the dependence on A 
and the right truncation point R. Note that the latter dependence is due to the maximum in Eq. ^ that is 
defined over the set St^, the set of all states that are reachable from a state in St a by at most R transitions. 
We have 



t+A 



A{t')dt'. 



(12) 



Here, we propose to first choose a desired right truncation point R* and then find a time-step A such 
that Eqs. ([7]) and ([H) hold. We perform an iteration where in each step we systematically choose different 
values for A and compare the associated right truncation point R with R* . Since /Xr*,a is monotone in 



A this can be done in a binary search fashion as described in Algorithm |l(a)| . We start with the two 
bounds A" = and A+ = ^^ax — The function FindMaxState(A,/?*) finds a state x"^'^ such that for all 
time-points t' € [t,t + A) we have 



£ ay(.^'^^^^') > max £ aj{x',t'). 

j=i j=i 

The choice of x"^^ also determines the uniformization rate 



(13) 



A(0 = I«y(^"^^0• 

7=1 

It immediately follows from Eq. (fT3] ) that Eq. ^ holds. In Section I4.2[ we discuss how to find A 
efficiently by selecting a state a:™^", while avoiding that the uniformization rates A{t') are chosen to be 
very large. 
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Input 



Output 



A, x" 



Global State space S, 



1 A+ := fmax — f ; //upper bound for A 

2 R:=0; 

3 x""^" :=FindMaxState(A+,^*); 

4 Mfi* .A+ : = ComputeParameter (f , f + A+ , jc™ 

5 ^+ := FoxGlynn(;Us._A+,£); 

6 if /?+ < ^* then 

7 A:=A+; 

8 else 

9 7?" := 0; A" := 0; //lower hound for A 
while ^ 7^ ^* 

A:= 



10 
11 
12 
13 
14 
15 
16 
17 
18 
19 



A+-A- 



Mr*, A := ComputeParameter (A, 7?*); 
R : = FoxGly nn (hr- ^a, £); 
if/? <R* <R 

R+ := R;A+ := A; 
e\seifR<R* <R+ 

R ■=R;A- := A; 
endif 
endwhile 



20 endif 



(a) The step size A is determined in a binary-search fash- 
ion. 



Input 



Output 



fo, fmax, p'°,e,R* 



p \ 



., p'' 



Global State space S, 



1 tcur '—to', 
3 tnpxt • — ^ru. 



: Algorithm 1{R*, tcur , f max , e ) 
f A; 



4 11 := ComputeParameter , f,ie.« , J^™^" ) ; 

5 while tcur < fmax 

6 / := 1 

7 while / < R* 

8 Compute v''' (jc); //DTMC probabilities 

9 Compute IPP N probabilities ; 
Accumulate p'""'{x); //CTMC probabilities 
i '.= ( + 1; 

endwhile 

tcur • — tficxt, 

(A,jc™^'') — Algorithm l(/?*,w,fmax,e); 

^maxv 



10 

11 

12 
13 
14 
15 

16 M • = ComputeParameter [tcur , tnext , 

17 endwhile 



(b) The complete algorithm. 



Fig. 1: Algorithms 

The function ComputeParameter + AjX™^") now computes the integral /i/j* a using a;"^^". If pos- 
sible we compute the integral analytically, otherwise we use a numerical integration technique. The 
function FoxGlynn(M,£) computes the right truncation point of a homogeneous Poisson process with 
rate pL for a given error bound e, i.e. the value R that is the smallest positive integer such that 

^ u' 

h " 

For the refinement of the bounds and A+ in hnes 13-17 we exploit that R is monotone in A. 



4.2 Determining the maximal rates 



The function FindMaxState(A,7?*) in Algorithm 1 1 (a) | finds a state x™^^ such that its exit-rate is greater 
or equal than the maximal exit-rate a()(x,?') = LyLi CCj{x,t') over all states x in St^R*. In principal it 
is enough to find a function A(?') with this property, for instance the function max;ci=5^^, ^7=i (^ji^jt')^ 
but this function may be hard to determine analytically and it is also not clear how to represent such a 
function practically in an implementation. Selecting a state x™^" and defining A(f') to be the exit-rate of 
this state solves these problems. 
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We now present two ways of implementing the function FindMaxState. 

a) For this approach we assume that all rate functions increase monotonically in the state variables. 
This is, for instance, always the case for models from chemical kinetics. We exploit that the change 
vectors are constant and define for each dimension k G {I, . . . ,n} 

je{l,...,m} 

where wjk is the ^-th entry of the change vector wj. For the set 5,.o we compute, the maximum 
value for each dimension k G {I, . . . ,n} 

j^^" := maxyk. 



We now find the state which is guaranteed to have a higher exit-rate than any state in 5, for 
all time-points in the interval [t,t + A) as follows, 

-^max ,,max , n* ...max 

--Jk +^ '^k ■ 

It is obvious that the state variables x™^" are upper bounds for the state variables appearing in 
St.R*. Then, since all rates increase monotonically in the state variables, we have that the exit- 
rate of x™^" = (x^"™, . . . ,xjf^'') must be an upper-bound for the exit-rates appearing in St.R* for all 
time-points. 

b) The first two moments of a Markov population model can be accurately approximated using the 
method of moments proposed by Engblom 16]. This approximation assumes that the expectations 
and the (co-)vaiiances change continuously and deterministically in time and it is accurate for 
most models with rate functions that are at most quadratic in the state variables. We approximate 
the means Ek{t') := E[Xk{t')] and the variances ol{t') := VAR[Xk{t')] for all k e {\,...,n}. For 
each k, we determine the time instant t £ [t,t + A) at which Ek{t) + £■ Ok{t) is maximal for some 
fixed i. We use this maximum to determine the spread of the distribution, i.e. we assume that 
the values of X{t') will stay below x™^" := Ek{t) + £■ Ok{i) with high probabihty. Note that a 
more detailed approach is to consider the multivariate normal distribution with mean E[X{t')] and 
covaiiance matrix COV[X {t')]. But since the spread of a multivariate normal distribution is difficult 
to derive in higher dimensions, we simply consider each dimension independently. We now have 
^max _ (^^max^^ ,x^^^). If during the analysis a state is found which exceeds x™^" in one dimension 
then we repeat our computation with a higher value for £. To make this approach efficient, i has 
to be chosen in an appropriate way. Our experimental results indicate that for two-dimensional 
systems the choice £ = 4 yields the best results. 



4.3 Complete algorithm 

Our complete algorithm now proceeds as follows (see Algorithm |l(b)| ). Given an initial distribution p^''' 
with finite support ^o^o, a time-bound t™'^^, thresholds 5 and e, and a desired right truncation point R*, 
we first set f := 0. Now we compute a time-step A and the state x""^" using Algorithm 1 1 (a) | with inputs 
R*, t, f™^", and e. We then approximate the transient distribution p'^^ using an on-the-fly version of the 
bounding approach 13, where the state space is dynamically maintained and states with probability less 
than 5 are discarded as described above. For the rate function A{t) we use the exit-rate of state x™^". 
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When computing DTMC probabilities, we use exact formulas for the first two terms IS of the sum in 
Eq. ^ and lower bounds, given by Eq. (ITOl) . for the rest. This gives us the approximation p'^^ with 
finite support St+Afi. We now set ? := ? + A and repeat the above step with initial distribution p' until we 
have t = 



5 Case Studies 

We implemented the approach outlined in Section|4]in C++ and ran experiments on a 2.4GHz Linux ma- 
chine with 4 GB of RAM. We consider a Markov population model that describes a network of chemical 
reactions. According to the theory of stochastic chemical kinetics [Sj, the form of the rate function of a 
reaction depends on how many molecules of each chemical species are needed for one instance of the 
reaction to occur. The relationship to the volume has been discussed in detail by Wolkenhauer et al. ifTTl . 
If no reactants are needec]!, that is, the reaction is of the form ©—)■... then aj{x,t) = kj • where 
kj is a positive constant and V{t) is the volume of the compartment in which the reactions take place. 
If one molecule is needed (case 5,- — )■ . . .) then aj{x,t) = kj -Xi where Xi is the number of molecules of 
type Si. Thus, in this case, aj{x,t) is independent of time. If two distinct molecules are needed (case 
Si + Se^ ■■ .) then aj{x,t) = A. . . 

All these theoretical considerations are based on the assumption that the chemical reactions are 
elementary, that is, they are not a combination of several reactions. Our example may contain non- 
elementaiy reactions and thus a realistic biological model may contain different volume dependencies. 
But since the focus of the paper is on the numerical algorithm, we do not aim for an accurate biological 
description here. 

We conduct experiments with two reaction networks. The first one is a simple gene expression 
(described in Ex. [T])- The second one is a gene regulatory network, called the exclusive switch |[T2l . It 
consists of two genes with a common promotor region. Each of the two gene products Pi and P2 inhibits 
the expression of the other product if a molecule is bound to the promotor region. More precisely, if 
the promotor region is free, molecules of both types Pi and P2 are produced. If a molecule of type Pi is 
bound to the promotor region, only molecules of type Pi aie produced. If a molecule of type P2 is bound 
to the promotor region, only molecules of type P2 are produced. No other configuration of the promotor 
region exists. The probability distribution of the exclusive switch is bistable which means that after a 
certain amount of time, the probability mass concentrates on two distinct regions in the state space. The 
system has five chemical species of which two have an infinite range, namely Pi and P2. We define the 
transition classes Xj = {Gj,Wj, aj), 7 G {1, . . . , 10} as follows. 

• For 7 G {1,2} we describe production of Pj by Gj = {x G | X3 > 0}, Wj = ej, and aj{x,t) = 
0.5 -xi,. Here, X3 denotes the number of unbound DNA molecules which is either zero or one and 
the vector ej is such that all its entries are zero except the j-th entry which is one. 

• We describe degradation of Py by Gy+2 = {x G | xy > 0}, wy+2 = —^j, and aj+zi^, t) = 0.005 -xy. 
Here, xy denotes the number of Py molecules. 

• We model the binding of Py to the promotor as Gy+4 = {x G | X3 > 0,xy > 0}, wy+4 = —ej — 
^3 + ey+3, and ay_|_4(x,f) = (0.1 — -t) -xy •X3 for t < 3600. Here, xy+3 is one if a molecule of 
type Py if bound to the promotor region and zero otherwise. 



Typically, reactions requiring no reactants are used in the case of open systems where it is assumed that the reaction is 
always possible at a constant rate and the reactant population is not explicitly modelled. 
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Case study 


FindMaxState 




Total 


Ex. time 






Poisson% 




R* 






mino/„ 


implementation 




error 










5 


4.69 • 10-"^ 


14 min 




95 


5 






10 


1.33 • 10^2 


10 min 




78 


22 




method a) 








33962 








15 


2.24 • 10^2 


5 min 




64 


36 


Gene 




20 


9.92- 10-2 


3 min 




41 


59 


expression 




5 


4.78-10-'* 


27 min 




95 


5 






10 


9.63 ■ 10 


14 mm 




77 


23 




method b) 








33130 








15 


4.08- 10-2 


10 min 




58 


42 






20 


7.73 - 10-2 


7 min 




41 


59 






5 


2.38 - 10-^ 


21 min 




80 


20 






10 


1.63 - 10-5 


29 min 




75 


25 




method a) 








1740 








15 


2.51 - 10-5 


68 min 




47 


53 


Exclusive 




20 


3.32-10-5 


2h 




38 


62 


switch 




5 


3.56-10-*^ 


17 h 




89 


11 




method b) 


10 


1.55-10-4 


3h 




78 


22 










1752 










15 


6.51 - 10-4 


1.5 h 




59 


41 






20 


1.71-10-3 


1 h 




42 


58 



Table 1 : Results of the analysis of case studies. 



• For unbinding of Pj we define Gy+e = {x G N | xy+3 > 0}, wj+e = ej + es — and ay+6(x,0 = 
0.005 ■Xj+3. 

• Finally, we have production of Pj if a molecule of type Pj is bound to the promotor, i.e., Gy+s = 
{x G N5 I A;y+3 > 0}, wy+8 = ^j, and aj+s{x,t) = 0.5 •xy+3. 

Note that only the rate functions and a(„ which denote the binding of a protein to the promotor 
region, are time-dependent. This is intuitively clear since if the cell volume grows it becomes less likely 
that a protein molecule is located close to the promotor region. We started the system at time f = in 
state (0,0, 1,0,0) with probability one and considered a time horizon of t = 3600. For the simple gene 
expression system (Example [B we started at time t = in state (0,0) and considered the same time 
horizon. Table [T] contains the results of our experiments. The first column refers to the system under 
study and the second one shows the variation used to implement the method FindMaxState which we 
suggest in Section 14.21 The third column lists the different values for right truncation point R*. We 
list the total error at time ?max in the fourth column (see Eq. (fTTI)). Program execution time is given in 
the fifth column and the sixth column with heading l^j contains the maximal size of the set St^R* that we 
considered during the analysis. The next two columns describe the percentage of the total probability loss 
due to the bounding approach (min%) and due to the truncation of the infinite sum in Eq. Q (Poisson%). 
The two percentages in one row do not sum up to one since we store only states that have significant 
probability (w.r.t threshold 5), which is the third error source. However, this lost portion is negligible for 
the two systems that we consider. For our implementation we kept the input e = 10-^" of Algorithm ] 1(a) 
fixed. 
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5.1 Discussion 

We now discuss the effect of the different input parameters on the performance of our algorithm and start 
with the implementation of the method to approximate For both systems the method "b" is less 

effective than method "a" (see Section l4!2l) . Method "b" gives larger uniformization rates than method 
"a", which leads to slower execution times. Notice that the execution time grows when we use method 
"a" for the exclusive switch system when we choose the larger values for R*. This is due to the fact that 
it always finds a state jc™^" without taking expectations and covariances into consideration. This results 
in large over-approximations for such a bi-stable system. The effect of the choice between methods "a" 
and "b" on the accuracy is not completely clear, both methods provide the same order of the probability 
loss for the simple gene expression system. For the second case study method "a" provides tighter error 
bounds for larger values of R*. 

In the Table[T]we show results obtained with 5 = 10^'^. Naturally, choosing a lower threshold results 
in larger execution times but one can gain a deeper exploration of the state space. This fact can also be 
used to obtain a coarse solution for certain system by setting 5 = 10^^, for instance. 

The effect of the choice of R* is most interesting. Choosing a larger value for R* means that more 
summands on the right-hand side of Eq. (flOl ) have to be approximated using the bounding approach. This 
decreases the accuracy of the algorithm since the larger time steps A are conducted and one obtain coarse 
approximation. However it reduces the running time since fmax can be covered using fewer iterations. 
Notice that the percentage of the probability loss due to truncation of the infinite sum in Eq. ^ grows 
when R* is chosen to be large. The reason is that we compute only first 3 exact terms in the sum and 
remaining terms are approximations. Thus the choice of R* determines the compromise between running 
time and accuracy. 

6 Conclusion 

We have presented an algorithm for the numerical approximation of transient distributions for infinite 
time-inhomogeneous Markov population models with unbounded rates. Our algorithm provides a strict 
lower bound for this transient distribution. There is a trade-off between the tightness of the bound and 
the performance of the algorithm, both in terms of computation time and required memory. 

As future work, we will investigate the relationship between the parameters of our approach (trun- 
cation point, the significance threshold 5, the method by which we determine the rate of the Poisson 
process), the accuracy and the running time of the algorithm more closely. For this we will consider 
Markov population models with different structures and dynamics. 

References 

[1] A. Arkin, J. Ross & H. H. McAdams (1998): Stochastic Kinetic Analysis of Developmental Pathway Bifur- 
cation in Phage X-Infected Escherichia coli Cells. Genetics 149, pp. 1633-1648. 

[2] M. Arns, P. Buchholz & A. Panchenko (2010); On the numerical Analysis of Inhomogeneous Continuous 
Time Markov Chains. INFORMS Journal on Computing 22, pp. 416-432, doif lO . 1287/i joc . 1090 . 
103571 

[3] G. Ciardo (1995): Discrete-time Markovian stochastic Petri nets. Kluwer. 

[4] F. Didier, T. A. Henzinger, M. Mateescu & V. Wolf (2009): Fast Adaptive Uniformization of the Chemical 
Master Equation. In: Proc. of HIBI pp. 118-127, doi ^lO . 1109/HiBi . 2009 . 23l 



A. Andreychenko, P. Crouzen, L. Mikeev, V. Wolf 



15 



[5] N.M. van Dijk (1992): Uniformization for nonhomogeneous Markov chains. Operations research letters 
12(5), pp. 283-291, doi jlO .1016/0167-63^7192) 9008'6^ : 

[6] S. Engblom (2006): Computing the moments of high dimensional solutions of the master equation. Appl. 
Math. Comput. 180, pp. 498-515, doi jlO . 1016/ j . amc . 2005 . 12 . 032] 

[7] D. T. Gillespie (1976): A General Method for Numerically Simulating the Time Evolution of Coupled Chem- 
ical Reactions. J. Comput. Phys. 22, pp. 403^34, doi:10Ti01670 02T£999TT7 6) 900 41-3[ 

[8] D. T. Gillespie (1977): Exact Stochastic Simulation of Coupled Chemical Reactions. J. Phys. Chem. 81(25), 
pp. 2340-2361, doi jlO .1021/ j 10054 0a008[ 

[9] W. K. Grassmann (1990): Computational methods in probability theory. In D. P. Heyman & M. J. Sobel, edi- 
tors: Stochastic Models, chapter 5. Handbooks in Operations Research and Management Science 2, Elsevier, 
pp. 199-254, doi:1 0.1016/S0 927-0507 (05) 80169-0. 

[10] A. Jensen (1953): Markoff chains as an aid in the study of Markoff processes. Skandinavisk Aktuarietidskrift 
36, pp. 87-91. 

[11] J. F. C. Kingman (1969): Markov Population Processes. Journal of Applied Probability 6(1), pp. 1-16. 

[12] A. Loinger, A. Lipshtat, N. Q. Balaban & O. Biham (2007): Stochastic simulations of genetic switch systems. 
Phys. Rev. E 75(2), p. 021904, doi ^lO . 1103/PhysRevE . 75 . 021904] 

[13] A. P. A. van Moorsel & K. Wolter (1998): Numerical Solution of Non-Homogeneous Markov Processes 
through Uniformization. In: Proc. of the European Simulation Multiconference - Simulation. SCS Europe, 
pp. 710-717. 

[14] E. de Souza e Silva & P. M. Ochoa (1992): State Space Exploration in Markov Models. In: SIGMETRICS. 
pp. 152-166, doi: 10 . 1145/133057 . 1331 00. 

[15] W. J. Stewart (1995): Introduction to the Numerical Solution of Markov Chains. Princeton University Press. 

[16] M. Thattai & A. van Oudenaarden (2001): Intrinsic noise in gene regulatory networks. PNAS, USA 98(15), 
pp. 8614-8619,doi ]10 .1073/pnas .151588598] 

[17] O. Wolkenhauer, M. Ullah, W. Kolch & K. Cho (2004): Modeling and Simulation of Intracellular Dynamics: 
Choosing an Appropriate Framework. IEEE Transactions on NanoBioscience 3(3), pp. 200-207, doi: 10 . 
1109/TNB.2004 . 833694] 



